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The asymmetrically forced, damped Duffing oscillator is intro- 
duced as a prototype model for analyzing the homoclinic tangle of 
symmetric dissipative systems with symmetry breaking disturbances. 
Even a slight fixed asymmetry in the perturbation may cause a sub- 
stantial change in the asymptotic behavior of the system, e.g. transi- 
tions from two sided to one sided strange attractors as the other pa- 
rameters are varied. Moreover, slight asymmetries may cause substan- 
tial asymmetries in the relative size of the basins of attraction of the 
unforced nearly symmetric attracting regions. These changes seems 
to be associated with homoclinic bifurcations. Numerical evidence in- 
dicates that strange attractors appear near curves corresponding to 
specific secondary homoclinic bifurcations. These curves are found 
using analytical perturbational tools. 
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Introduction 



The 



forced and damped Duffing oscillator: 



x -\-sS x —x + x 3 = £7cos(cut), (x,t) G 3? 1 x 
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has served as a prototype model for investigating low dimensional chaotic 
behavior in enumerate publications (see [16, 22, 37] and references therein). 
Its significance lies in its simple "typical form" which appears in many ap- 
plications. Indeed, the unperturbed Duffing oscillator represents the normal 
form for Hamiltonian systems with Z 2 symmetry [13]. Thus, whether dif- 
ferent types of perturbations lead to substantially different dynamics is of 
mathematical and physical significance. The perturbation of (1) has two 
specific properties - it has no non-linear terms in i,i and it is symmetric; 
(1) is invariant under x — > —x } t — > t + tt/lo. 

Numerical simulations suggest that the inclusion of nonlinear dissipation 
term in the perturbation does not alter the qualitative behavior of the forced 
system [29]. Namely, no new bifurcation sequences or new types of attractors 
appear, though the location of the various bifurcation curves of (1) changes. 
The effect of asymmetric potentials has been investigated when the forcing 
is adiabatic, see [25, 3] and references therein. In this paper, we examine 
the effect of asymmetric forcing on the Duffing oscillator by introducing the 
asymmetrically forced, damped, Duffing oscillator (AFDO): 

x +6 x -x + x 3 = (x - /3x 2 )7cos(cut), (x,t) G 3? 1 x (2) 

which contains the asymmetry perturbation parameter f3. Here we show that 
the inclusion of the physically typical asymmetric forcing perturbations alters 
the qualitative behavior of the system in some range of parameter values. 

Theoretical and numerical investigations of forced and damped systems 
with homoclinic tangle is problematic since these may attain Strange Attrac- 
tors (SA) 1 and periodic sinks simultaneously. Moreover, the existence of the 
SA is extremely sensitive to changes in parameter values. The existence of 
Newhouse sinks near homoclinic tangencies implies that small changes in the 
parameters may destroy the SA. These observations are reflected in the dif- 
ficulties of proving the existence of SA in such systems (see reviews [16, 37]). 
Analytical results regarding the existence of SA [1], their basins of attraction 
[5], and the construction of unique natural invariant measure (the SRB mea- 
sure) [2], has been recently published for the Henon map. These proofs are 
in the strong dissipation limit, for which the strange attractor appears as a 
one dimensional attractor multiplied by a cantor set. Some of these results 
may be applied locally to neighborhoods of homoclinic tangencies [23]. 

1 These are attractors with sensitive dependence on initial conditions, i.e. attractors 
which have a dense orbit with positive Lyapunov exponent (see [9]). 
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It follows that SA are expected to emerge near homoclinic tangencies. 
We use analytical tools for locating primary 2 homoclinic tangencies (the Mel- 
nikov analysis, [16]) and secondary 3 homoclinic tangencies (the SMF [33]). 
See [34, 17, 18, 4] for other works on the subject of multi-pulse homoclinic or- 
bits in other settings. This presents the first application of the SMF to a dis- 
sipative system. Since homoclinic bifurcations are considered an important 
source of structural instabilities of dynamical systems [27, 11, 12, 14, 19, 26], 
their location in parameter space should indicate regions in which dramatic 
structural changes appear. Clearly higher order tangencies exist as well, and 
finding them all is a useless mission, in particular in view of Newhouse work. 
The philosophy here is that not all homoclinic tangencies have the same sig- 
nificance: primary tangencies are more important then secondary, secondary 
more then third order etc.. Thus there is a sense in locating the bifurcation 
curves of the lower order homoclinic tangencies. This approach is backed 
up by the TAM (topological approximation method [31, 32]), which asserts 
that many features of the dynamical system are determined already by the 
characteristics of the primary and secondary homoclinic orbits. The TAM 
was developed for non-dissipative systems and has been recently generalized 
to dissipative systems [20]. 

The last part of this work consists of a numerical search for SA at param- 
eter values which are close to the analytically predicted bifurcation curves. 
SA have been observed in various systems exhibiting homoclinic chaos, in- 
cluding the forced and damped Duffing oscillator [16], the Henon map [19], 
and the forced and damped cubic potential [21]. In the latter work the 
correspondence between the appearance of homoclinic tangencies of specific 
character and SA has been noted, correspondence which seems to persist for 
the AFDO. 

This paper is ordered as follows: In section 2 we present the basic phase 
space structure of the AFDO, the Melnikov analysis and the bifurcation 
curves for primary and secondary homoclinic bifurcations. Numerical ev- 
idence suggesting the existence of SA near specific homoclinic bifurcation 
curves is presented in section 3, the typical size and shape of the basin 

of attractions of the attractors. Conclusions and a discussion are presented 
in section 4. 

2 These are one-loop homoclinic orbits which are O(s) close to the unperturbed homo- 
clinic orbits for t £ (— oo, oo). 

3 These are two-loops homoclinic orbits which are O(s) close to the unperturbed homo- 
clinic orbits for t £ (— oo,t ], [ti, oo), see also section 2. 
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2 Templates of the homoclinic tangle 



2.1 Basic properties of the AFDO 

Introducing the phase space coordinates (x,y) £ 3? 2 , one rewrites equation 
(2), as: 

x = y, (3) 
y = x — x 3 ± (x — f3x 2 )e~{ cos(cut) — s8y , 

Physically, 8 represents the dissipation (the damping), 7 the amplitude of 
the forcing, u the frequency and f3 the asymmetry disturbances. These pa- 
rameters are real and by symmetry may be taken to be non negative, e is 
a "perturbation scaling parameter", assumed to be small. For £ / 0, there 
are two differences between the AFDO (2) and the Duffing oscillator ( I). 
The substantial difference is that (2) includes the asymmetry parameter, f3. 
The second difference is that for convenience, with no loss of generality, the 
symmetry x — > —x } t — > t ± tt /u of (I) is replaced by the symmetry x — > — x 
for f3 = in (2), hence the origin is fixed for all e,/3. f3 ^ corresponds to 
symmetry breaking disturbances. 

The unperturbed system corresponds to the integrable Hamiltonian sys- 
tem with a symmetric quartic potential: 

? 4 

V(x) = - X -+ X -, (4) 
and with the Hamiltonian function (energy): 

9 9 9 4 

<i t it r f* r f* 

H(x,y) = y Y + V(x) = y Y - Y + T . (5) 

The unperturbed system, which is identical to that of the unperturbed Duff- 
ing oscillator, has three equilibrium points: two centers at (x,y) = (±1,0), 
and a saddle at (x, y) = (0, 0). The saddle point is connected to itself by two 
homoclinic orbits, with periodic orbits nested within and around them. The 
period of the unperturbed periodic orbits, P(H) } has the following asymp- 
totic expansion near H = (exact formulae for all H are available [16, 37]): 

ln(i|)(l±0(i7)), H^O- 
P(H) = I (6) 
21n(^)(l±0(#)), #-0 + . 
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In the unperturbed system the stable and the unstable manifolds of the 
saddle point (0, 0) coincide. For 8 > 0, and 7 = 0, the unstable manifold of 
the saddle point near the origin falls into the two sinks created near (±1, 0). 
As for the Duffing oscillator, it may be proved that for sufficiently small 
values of 7 the closure of the unstable manifold (which contains the saddle 
and the sinks) is an attracting set of (3). 

A Poincare map in time is used to simplify the phase space portrait for 
the time dependent system (7 ^ 0). Keeping 8 > 0, and increasing 7, the 
following scenario occurs on both sides of the fixed point; for small values 
of 7, the Poincare map is topologically equivalent to the Poincare map with 
7 = 0, which is structurally stable. As 7 increases, resonance bands of 
higher period and higher amplitudes are created. As 7 is further increased, 
in addition to the resonances, a homoclinic bifurcation occurs, after which 
the stable and the unstable manifolds of the saddle point of the Poincare map 
intersect in transversal homoclinic orbits. The presence of these orbits implies 
the existence of a complicated non-wandering Cantor set which possesses 
infinitely many unstable periodic orbits of arbitrary long period as well as 
bounded non-periodic motions. The Smale - Birkhoff Homoclinic Theorem 
implies that in this case the system has chaotic dynamics. When f3 = 0, 
as in the forced Duffing oscillator (1), the sequence of bifurcations described 
above occurs simultaneously on both sides of the fixed point. When f3 > 
this changes as described below. 




o) m b 

a) F 1 (u), F 2 (u) } b) The ratio 

Figure 1: The Melnikov function amplitude. 
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The Melnikov function, M(t + 6/ 00), measures the signed distance be- 
tween the stable and the unstable manifolds of a hyperbolic fixed point (up to 
a multiplication by a constant). This distance is measured at the Poincare 
section cot = 6 } and to represents a parameterization along the unstable 
manifold. For the AFDO (equation (3)) two Melnikov functions are defined: 
M r (t ; 7,cu, /3, 8) = M r (t ) (respectively M/(t )) measures the signed distance 
between the right (respectively left) branches of the stable and unstable man- 
ifolds. These functions are given by: 



M rtl (t ;i,u,P,6) = / (yx(l - /3x)-fcos(iot) - 6y ) \ {q o (t)rUt+to) dt 

J —GO 

4 

= 'ysm(ut )F rt i(u,(3) - -8, (7) 

where q°(t) r j are the right and left unperturbed homoclinic orbits of the 
system: 

q°(t) r ,i = ±(v^secht, -v^sechttanh t), (8) 



an> 



d 



F r {u,P) = [FiH-^H], (9) 

F(u^) = [^H + ^H], (10) 



Fl(io) = 7TC0 2 CSch(^), (11) 



F 2 (u) = — ttuj(1 +cu 2 )sech(— ). (12) 



a/2 2 ,TTLO, 

— iruj(l + u ) sech( — 
3 2 

Figure 1 shows Fi(u), F 2 {to)^ and the relation ^ . Notice that Fi(uj) and 
F 2 (io) are non-negative for all u. 



2.2 Primary homoclinic intersection points 

For e sufficiently small simple (respectively degenerate) zeros of the Melnikov 
function imply primary homoclinic transverse intersections (respectively tan- 
gencies) of the stable and the unstable manifolds of the hyperbolic fixed point 
[16]. Requiring M(t ) = M'(t ) = 0, it follows from (7) that primary homo- 
clinic bifurcations occur near 



e-f r: i(8;u,/3) = e8 



3[F 1 (uj)t/3F 2 (uj)] 



e8B±(u,,/3). (13) 
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3 7 IFxH - (3F 2 {lo)\ ^[F^lo) + I3F 2 {lo)} 
< o < 



(14) 



4 4 
then the left branches of the stable and the unstable manifolds intersect, 
while the right branches do not, see for example figure 2 (the intersections 
on the right hand side of this figure correspond to secondary homoclinic 
points, as described in 2.3). Similarly, if 



S < 



^IF^uj)- f3F 2 (uj)\ 



(15) 



then the stable and unstable manifolds intersect on both left and right sides 
of the hyperbolic fixed point, see for example figure 7. It follows that the 
parameter space is divided to three regions (see figures 3, 4): 

I For j < Rq(lo } /3) there are no primary intersections of the stable and 
the unstable manifolds. 

II For Rq(u,/3) < j < Rq(u,/3) the stable and the unstable manifolds 
have primary intersection points on the left side of the saddle point, and 
do not have primary intersection points on the right side (asymmetric 
behavior). 
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Ill For j > Rq(u,/3) primary intersections of the stable and the unstable 
manifolds occur both on the left and the right sides of the saddle point 
(like in (1), but in an asymmetric manner for (3^0). 



Si 



mce 



i?o>,/3) 



1 



(16) 



the relative size of region II depends on the values of x = f3 , and may 
be derived from the graph of r(x) (hgure 5) for the corresponding f3 and u 
values; Fixing f3 ^ 0, the relative size of region II varies with u as described 
below. First, notice that is bounded from below by a positive constant 

c ~ 0.7, and that it grows monotonically (in fact asymptotically linearly) 
with to (see hgure 1). Therefore, and since r(x) = 1 at x = 0, oo only, it 
follows that for any finite non- vanishing u value region II is of non- vanishing 
measure (see for example figures 5 and 4). Since x grows monotonically with 
> /3c , the relative size of region II increases with u up to the 



cu, and x 



8 



S/y 



Figure 4: Primary homoclinic bifurcation curves (in space). 

threshold value io*(/3) at which = ^ (i-e. x({3,uj*({3)) = l) 4 . For x > 1, 

r(x) strictly increases with x, thus for > io*(/3) the relative size of region 
II decreases as u increases. At x = 1, r(x) = 0, namely Rq(uj,/3) — > oo as 
— )■ cu*. For these values, regions I and II occupy most of the parameter 
space, and region III shrinks till it disappears, to order e, at u = u*(f3). 




U I <L J t U U Lf.J I I .J d 

a b 

a) ^(u,f3 = 0.1), see eq. (13 ), b) r(x), see eq. (16). 

Figure 5: The relative size of region II. 
4 If j3 > niai^( ^H ) = then does not exist, and may be considered as infinite. 
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2.3 Secondary homoclinic intersection points. 

First, we describe geometrically what are the secondary homoclinic intersec- 
tion points, what are their transition numbers and what are the structural 
indices of a homoclinic tangle. Then, we present the analytical (perturba- 
tional) method for finding these structural indices. 

Consider values of the dissipation parameter 8 for which the Melnikov 
function M/(t ; //, 8) (respectively M r (t ; //, 8)) has two simple zeros (see sec- 
tion 2.2). Denote the corresponding PIPs (Primary homoclinic Intersection 
Points), ordered by the direction of the unstable manifold, by pl (respec- 
tively pr ) and ql (respectively qr ) see figure 6. Also, denote their or- 
dered images under the Poincare map F by pk^qk (respectively pr 8 ,gr 8 ), 
i = 0, ±1, ±2, ±3, i.e. F l (p! ) = pli and so on. The areas enclosed by 
the segments of the stable and the unstable manifolds connecting two suc- 
cessive PIPs are called lobes. Denote the lobes enclosed by segments of the 
stable and the unstable manifolds connecting pk^qk by Dl{ } and the ones 
that between qli 7 pli + i } by £7;, when again, Dli = F l (Dlo), Eli = F l (Elo) 
(the equivalent notation is used for the right side), see figure 6. 
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Figure 6: Illustration of Lobes, PIPs and SIPs. 
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If Elj H Dl ^ (respectively Erj fl Dr ^ 0) for some non negative 
integer j, or if Dlk+i H Er ^ (respectively Dr^+i fl El ^ 0) for some 
non negative integer k } then there exist Secondary Intersection Points (SIPs) 
in these intersections. The integers j, k are the transition numbers of the 
corresponding SIPs [8]. The minimal transition numbers (the minimal inte- 
gers j, k) for which this happens on the left side (respectively right side) of 
the hyperbolic fixed point are called the the structural indices £u,£i r (respec- 
tively £ rr} £ri) of the homoclinic tangle [8]. Namely, these structural indices 
are exactly the transition numbers of the secondary homoclinic points which 
belong to the first intersection of the corresponding lobes. For example, in 
figure 6 the structural indices are: £\\ = 1, £\ r = 1, £ r \ = 1 and £ rr = 2. 
Each such structural index imposes minimal complexity for the structure of 
the homoclinic tangle. For example, the length growth rate of line segments 
along the unstable manifold (the topological entropy) increases as the indices 
decrease. Thus, dividing the parameter space to regions according to the val- 
ues of these indices corresponds, approximately, to a plot of "level sets" of the 
topological entropy. The dividing curves correspond to secondary homoclinic 
bifurcations. Large SA which are not associated with primary homoclinic 
tangencies seems to appear only in regions in which at least one structural 
index is less or equal to 1 (see also [21]). 

The perturbational method for calculating the secondary homoclinic bi- 
furcation curves is described below. For simplicity it is presented specifically 
for the AFDO model. More generally, it may be applied to nearly Hamilto- 
nian dissipative systems, which satisfy some generic assumptions (see [20]). 

Consider the Secondary Melnikov Function - SMF (see [31, 32, 33] and 



P(H) is the period of the unperturbed periodic orbit with energy H } and 
H = on the separatrix. For sufficiently small e, simple zeros (respectively 
degenerate zeros) of (17) imply transverse secondary homoclinic intersections 




(17) 






< t < T, 



(19) 
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where [x] is the integer part of x, T = ^ is the period of the perturbation and 
s(t ) is either or 1, depending on the interval to which t belongs (equation 
(24) below). The structural index £ cd (c,d £ {h r }) i s defined to be the 
minimal transition number j c d(to } £)- For sufficiently small e this analytical 
definition of the structural index meets the geometrical definition described 
above [33]. It follows that typically 5 a change in the structural index may be 
found at a bifurcation point for equation (17). Indeed, under some generic 
conditions on h^{to,e\ [J-,8), the structural indices satisfy [33]: 

£ c d = jcd{tocd,£cd); c,d£{l,r}, (20) 

where (to cd} e cd ) (^cd small) are the solutions to the standard equations for a 
bifurcation point: 

hf{t Q ,e) = 0; 0<t o <T (21) 

= °- (22) 

defined in the appropriate time-interval for t\ : 

[£ c d + s(t 0cd )]T < t lcd (t 0cd ,e cd ) < [£ c d + s(t 0cd ) + 1]T; (23) 



s(t , 



cd 




'IT 



(24) 



where ti cd (t 0} e) is defined by equation (18). 

Typically, for e cd sufficiently small, one finds a sequence of two bifurcation 
values, e\ d < e 2 cd . The corresponding solutions (t l 0cd , e l cd ) ; i = 1,2 of (21) and 
(22) divide the parameter space to three regions: below the hyper-surface 
e = £^(w, 7, /3, <5, £ cd ) there are no SIPs, between the hyper-surfaces e = 
ej^w, 7, /3, <5, £ cd ) and e = £^(w, 7, /3, 6, £ cd ) two SIPs occur (see for example 
the intersections denoted by an arrow in figure 7), and above the hyper- 
surface e = £^(w, 7, /3, <5, £ cd ) two additional SIPs occur (see for example the 
intersections of the lobes above the origin in figure 7). 

Moreover, equations (17) - (22) may be brought to a simple form, as shown 
in appendix A. Using the asymptotic expansion for the period function P(H) 
(equation (6)), these equations may be solved analytically if both (3 = and 



5 Another source for changes in the structural indices are points of discontinuity of 
*i(*o), see [20]. 
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8 = 0, to find approximations to the secondary homoclinic bifurcation points 
(4,e' c d(*o))> c -> d e i l -> r }-> z e {^2}. The solutions for f3 = 0, 8 = may 
be used to solve these equations for f3 ^ or/and <5 ^ small by the use 
of asymptotic expansions in powers of f3 and 8 for i = 1, and in powers 
of /3~ and <5~ for z = 2. To get more accurate results, using P(H) instead 
of its approximate form, we use a Newton method, combined with a linear 
prolongation scheme (see appendix A and [20] for details). Plotting the 
bifurcation values, e 8 C(i (cu, 7, /3, 8, i c d)\ c, d £ 0,^}, i £ { 1 , 2} , obtained from 
the solutions to these equations, for 7, /3, 6 fixed and varying cu, gives the 
secondary homoclinic bifurcation curves in parameter space (cu,e), labeled 
by the structural indices, = 0, 1, 2, m < 00, as in figure 7a. 

A simple lower bound to the homoclinic bifurcation curves (compare with 
(30)) is given by: 



r p 



-1 / 2vr 



(4i + 1)) 



e cd (u,-f,/3,8,£ cd ) = < 
Therefore, using (6,7), we find 

e cd 2 (^7?/M,4i) > 



max to |M c (t )| 

P- 1 (^(4, + l)) i 

max to |M c (t )| 



c = d 



c ^ d 



(25) 



16e- 2 7(^+ 1 ) 



{| 7 F C -^|,|7^ + ^|}- 



max 



(26) 



These give a simple lower bounds on the secondary homoclinic bifurcation 
curves 6 . Moreover, geometrically these lower bounds correspond to the values 
of e for which the lobes may get involved in a 1 : (l c d + 1) resonance (see 
below). 



2.3.1 Comparison between numerical and analytical results 

The analytical method described above for finding the secondary homoclinic 
bifurcations is of a perturbational nature. Thus, as proved in [33], in the 

6 Note that since the approximation to leading order in H for the period function, P(H), 
is used here to calculate e c d, we get that e c d = s cc (see the appendix for more details). 
Hence, the curve s cc (£ cc = n) serves as a lower bound to all the eight secondary homoclinic 
bifurcation curves, related to the structural index £ cc = n , £^(£ C c = n ), with c, d £ {/, r} 
and n > 0. 
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Figure 7: Secondary homoclinic bifurcation curves and the homoclinic tangle. 



limit of small e values it is guaranteed to supply a good approximation to 
the actual bifurcation value. Here we examine how good of an approximation 
the analytical formulae supply for finite e values. Indeed, excellent agreement 
is achieved between the analytical predictions for the occurrence of SIPs and 
the numerical calculations of the stable and unstable manifolds for e as large 
as 0.3 and > 1; c, d £ 0,^}, see for example figure 7. In this figure the 
'x' at the (u:,e) parameter space indicates the parameter values for which 
the manifolds, presented on the right figure, are calculated. On the right 
figure, the corresponding near-tangency of the manifolds is indicated by an 
arrow. In fact, the larger the £'s, the larger the e values for which the zeroth- 
order approximation is found to be adequate. For example, for £ = 2 we find 
excellent agreement up to e 1. This is not surprising since large €s (and 
finite e) correspond to large cu's for which the Melnikov function coefficient 
becomes exponentially small, thus the effective perturbation is small. 

For £ cc = the agreement between the numerical and the analytical re- 
sults is not as favorable (notice that this is a finite e effect: letting e — > 0, 
with all other parameters held fixed, necessarily implies that £ — > oo); This 
is due to the passage of the manifolds through a 1 : 1 resonance relation 
between the periodic orbits inside the homoclinic loop and the forcing pe- 
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riod 2tt/uj; Namely, the energy level to which the manifolds are pushed 
by the Melnikov function is near the energy level for which a 1 : 1 reso- 
nance occurs 7 Now, the construction of the SMF uses the Whisker map (see 
[6, 8, 10, 30, 31, 32, 33, 40]) in which the motion of the interior orbits is ap- 
proximated by unperturbed periodic motion. This approximation fails near 
a 1 : 1 resonance. Hence, the analytical approximation for the £ cc = bi- 
furcation curve are inaccurate even for small values of e. Notice that in this 
limit to is varied with e, hence this observation is not contradictory to the 
SMF theorems which hold in the limit e — > with all other parameters held 
fixed [33]. Indeed, to avoid passage of the manifolds through a 1 : 1 reso- 



le > 1. For the outer indices £ c d,c ^ d the problem of 1 : 1 resonance was 
not encountered. 

3 Strange attractors. 

In this section numerical evidence for the existence of SA, and observations 
regarding their location in parameter space and their structural properties 
in phase space are presented. 

3.1 Numerical scheme for detecting SA. 

Simple numerical experiments showing Poincare maps of the AFDO (using 
DSTOOLS [7]) suggested that SA appear in the area of the parameter space 
related to the structural indices = 0; c, d £ {h r }- To investigate this 
subject more thoroughly the Lyapunov exponents of orbits of (3) were com- 
puted (see [38]). Viewing (3) as an autonomous system, each orbit has three 
Lyapunov exponents, one zero, one negative and the third may be either posi- 
tive or negative. A positive third Lyapunov exponent indicates the existence 
of a strange attractor (see [9]), while a negative third Lyapunov exponent 
indicates that the orbit is attracted to a periodic sink. 

An efficient stopping criteria for the Lyapunov exponents calculation is 
developed, using the distinction between SA or sinks with long transients 

7 Indeed, the 1 : m resonance relation for the periodic orbits of equation (3) is given by: 
P(H) = . Since by definition of t lcc (see equation (18)), P(s cc M c (t 0cc )) = ti cc -t 0cc , 
and by condition (23), t\ cc — to cc £ [^cc~, (^cc + 1)77); the manifolds "pass" through the 
1 : 1 resonance zone for l cc = 0. 



nance, e should satisfy the condition: e cc < 




This condition holds for 
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and simple sinks. First, to remove transient behavior, N{ n iterations ol the 
Poincare map 8 are calculated. Then, every Na iterations (Nn = 100) a line 
is fitted to the logarithm of the modulus of the last Na Lyapunov exponent 
values. The program stops if one of the following events occurs: 

1. The largest Lyapunov exponent is negative. Then there exist a periodic 
sink. The exact value of the negative Lyapunov exponent is not sought. 

2. The largest Lyapunov exponent is positive and the slope of the fitted 
line is nearly zero (up to an error of le — 6). 

3. The total number of Poincare map iterations exceeds 10000. In this 
case no decision is made regarding the existence or non-existence of an 
attractor. In practice the stopping criteria 1 and 2 occur before 10000 
iterates are computed. 

In case 2, when a positive Lyapunov exponent is detected, all the Lya- 
punov exponents are calculated. Hence, the Lyapunov dimension [24] of the 
strange attractor may be calculated. Note that the Lyapunov dimension, Dl 
, is an upper limit for the Capacity (or Box-Counting) dimension, D [24]. 

In the numerical experiments, 7, /3, 8 are fixed and e and u are varied along 
and near the secondary homoclinic bifurcation curves £' C(j (a;, 7, /3, 8, £ c d)i i = 
1,2; c, d £ {/, r}, of section 2.3. 

3.2 Windows of SA. 

For various parameter values, numerical evidence suggests the existence of 
SA in "windows" in the parameter space. These windows are aligned near the 
secondary homoclinic bifurcation curves which are related to the structural 
indices £u,£i r = 0, see figure 9 8, and the magnification of the windows in 
figure 9. While theoretically such regions should appear near all tangent 
bifurcation with arbitrary £ } we did not detect in our numerical search any 
SA near the bifurcation curves with £ > 0. This suggests that the size of the 
parameter regions for which SA appear decreases dramatically with £. 

8 Ni„ = 200 was found sufficient. 

9 Notice that in figure 8a (3 = hence e l u = e l rr and e\ r = e l rl) whereas in figure 8b 
£ \i i 1 £ rr an d e \ r 7^ £ w ' ^ ut > f° r clarity, s l rr , s l rl are not plotted. In figure 8c e l rr , s l rl are 
not defined for the specified u> values (since the Melnikov function M r (to) has no zeros). 
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Figure 9: Secondary homoclinic bifurcation curves and SA for the AFDO. 
Magnification of certain regions of figures a-c. 



18 



In fact, our perturbational methods for detecting the homoclinic tan- 
gencies associated with the structural indices In = are inaccurate, see 
section 2.3. Nonetheless, the predicted analytical bifurcation curves for 
l c d = 0; c, d £ still lie in the area of parameter space near which 

the actual bifurcation curves exist. Moreover, observe that the SA appear 
only in the region of the parameter space (cu,e), which is above the second 
secondary bifurcation curve sfi(uj, 7, /3, <5, In = 1) (see section 2.3). For e < 1, 
our predictions for this curve are accurate, hence, this curve may be consid- 
ered as a lower bound to the region in parameter space in which SA appear, 
see figures 8. 

Another feature of the SA windows is that they all seem to appear above 
a threshold value e(cu, 7, /3, 6) > 0.2. Namely, they do not seem to extend 
to the small e values to which some of the I = bifurcation curves extend. 
Numerical calculations of the stable and the unstable manifolds of the origin 
for the minimal e values for which strange attractors are found suggest that 
this curve is a specific homoclinic bifurcation curve; above this curve the lobe 
El\ intersects the lobe Dl at five or more (six, seven or eight) homoclinic 
intersection points. See for example figure 10. Notice that the homoclinic 
bifurcation curve £^(u;, 7, /3, <5, In = 1) in the (00, e) parameter space gives a 
lower bound to this curve. 




Figure 10: The stable and the unstable manifolds corresponding to SA. 

The structure of the SA that are obtained vary with the parameters. The 
main forms of the attractors which were found are described next. 
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3.3 The phase space structure of the SA. 

The observed SA have the following two distinct properties: 

1. The attractors may be one sided (i.e. the attractor is contained in 
the right/left half phase space plane) or two sided (with one or two 
components). 

2. The attractors may have strong dissipative features or may have nearly 
conservative features. 

The first property depends on the location of the parameter values with 
respect to the division to regions I, II and III. One sided SA may appear 
in region II near secondary homoclinic bifurcations or near the boarders 
between the regions, namely near primary homoclinic tangencies, see figure 
11. There, the transition between two-sided SA, denoted by 'T', and one 
sided SA, denoted by '0', is shown 10 . In fact, near the boarder between 
regions II and III three different SA may appear: one sided SA on the right 
half plane, two one-sided SA or a two-sided SA. Near the boarder between 
regions I and II left sided SA, coexisting with a sink on the right half plane, 
were observed. 

The second property seems to depend mainly on the ratio 6/7 and is 
roughly independent of the other parameters (in regions where SA exist); 
This is somewhat surprising since the area contraction per Poincare map is 
given by exp (— s62tt /u) - thus strong dissipation may be achieved for fixed 
(5/7 by increase of e, without essential changes in the structure of the SA. For 
^ < 1, the two sided strange attractors seems to have nearly conservative 
features of a chaotic region, see figure 12. These features persist in the 
window shown in figure 9a, and even when e = 1, though the attractor is 
more structured, it has "fat" regions in which no filamentation is observed. 
The positive Lyapunov exponent corresponding to this figure is log 2 ( 
0.1987, and the corresponding Lyapunov dimension is Dl ~ 1.9036. The SA 
indicated in figure 8a are of such structure. Most of these attractors have a 
positive Lyapunov exponent of about 0.2. The maximal observed deviation 
from this value is 0.02. Nearly conservative one-sided SA were not observed. 



10 This transition is not continuous in -; Between the value of - for which a two sided 

7 ' 7 

SA appears, and the value of ^ for which a one sided SA appears, there may be some 
values of — for which no SA appear. 
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Figure 12: Nearly-conservative strange attractor. 

When ^ is not small (6/7 = 0.2 is already in this region), the structure of 
the SA is associated with the folding of the unstable manifold as in strongly 
dissipative systems, see figure 13. The corresponding positive Lyapunov ex- 
ponent is log 2 ( 0.1741, and the Lyapunov dimension is Dl ~ 1.2975. 
Comparing with the above results for the nearly-conservative attractors, we 
observe that the Lyapunov exponent is less sensitive to the attractor's struc- 
ture than the Lyapunov dimension. The SA indicated in figure 9b are of 
strongly dissipative nature; some are two-sided with a positive Lyapunov 
exponent very close to 0.17, with maximal deviation of 0.03, and some are 
one-sided. In figure 14 such one-sided SA is shown, its positive Lyapunov 
exponent is log 2 ( 0.0862, and its Lyapunov dimension is Dl ~ 1.1453. 
The SA presented in figure 8c are all one-sided strongly dissipative (Henon- 
like) SA. r The values of the positive Lyapunov exponent are 0.08 ± 0.04, 
about half of the Lyapunov exponents of the two-sided SA. 

The Lyapunov exponent and dimension of the attractors seems to be 
quite robust. The dependence of the Lyapunov exponents on the values of e 
along a secondary homoclinic bifurcation curve related to a structural index 
^cd = is shown in figure 15. An example for the dependence of the Lyapunov 
dimensions on the values of e along such bifurcation curve is shown in figure 
16. The plunges in the figure correspond to parameter values for which no 
SA exist. 
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Figure 13: Two sided strange attractor of the AFDO, with structure referring 
to large values of the relation ^. 



In [20] symbolic dynamics of segments of the unstable manifold is con- 
structed for the AFDO, and for general dissipative systems which unfold 
homoclinic tangencies. From this symbolic dynamics a transfer matrix may 
be constructed for each set of the structural indices i c ^\ c,d £ {h r }- It fol- 
lows from [28] and [32], that log(A), where A is the modulus of the largest 
eigenvalue of the transfer matrix, gives a lower bound on the topological 
entropy of the Poincare map. The lower bound on the topological entropy 
for the AFDO with = 0; c, d £ 0,^}, corresponding to the region in 
parameter space where two sided SA appear, is log 2 (3.9231) = 1.9720. The 
lower bound on the topological entropy for the AFDO, corresponding to the 
existence of one sided SA is log 2 (3.6709) = 1.8761. 

The above results are consistent with the inequalities describing the rela- 
tions between topological entropy, entropy (Kolmogorov-Sinai invariant) and 
positive Lyapunov exponents (see [9] and [39]): 

h{p) < h top 

Hp) < E A « 

A 8 >0 

Where, p is an ergodic measure with compact support, with respect to a 
diffeomorphic map F, and A; are the positive Lyapunov exponents corre- 
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Figure 14: One sided strange attractor. 

sponding to a dense orbit of F . i.e. for the structural indices = we 
get: 

h(p) < Ai log 2 («i) < log 2 (A) < h top . 

Notice that the lower bounds obtained for the topological entropy are by 
order of magnitude larger then the corresponding Lyapunov exponents and 
that there is a very little difference between the one-sided and the two-sided 
cases. See [31] and [20] for the construction of the symbolic dynamics of the 
lobes, and details on how the transfer matrices and the lower bounds on the 
topological entropy may be calculated. 

3.4 Basins of attraction. 

Initial conditions may be attracted to the various attractors which exist in the 
phase space. Here, we distinguish between three types of attractors: those 
located entirely on the left (respectively right) half plane and those which 
are located on both sides of the y axis. We do not distinguish here between 
the basins of attraction of different sinks or SA, see [26, 36, 35, 15, 21] for 
detailed study of these issues. 

The flux of phase space area into the left /right sides may be calculated 
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Figure 15: The positive Lyapunov exponent variation. 

a) Samples ol e values along the bifurcation curve e = sj r (uj, 7, /3, 6; £i r ) with 
fj = 0, i = 0.05 and 4. =0. 

b) Zoom in on e values in the interval [0.2,0.25]. 

to hrst order in e by integrating the Melnikov function: 



l.r 



Mi, r (s)ds + / Mi, r (s)d. 



Sl 



" 3a; ' S ° 



si, M,, r (s) < Vs £ [0,T) 



1 . , 4<5 
s = — arcsml — 



/ 45 



7T 

),Sl = — 

to 



|^(2arcsin( T ^ 

■50, 



(27) 



+ 7r); s ^ 5 i- 



where determined by M^ r (s) < for s < s < Si £ [0, T\. Thus, on 

exp(2TveSn/u) , is swept into 



the n th iterate the initial phase space area, A 
the left/right side. In region I, M Ur {s) < Vs £ [0,T), and |A™| = |A^| 
to order e. Since F r < Fi it can be easily shown that |A™| > |A™| in regions 
II and III. Thus, for parameter values corresponding to these regions, the 
influx to the right side is always larger than the influx to the left side (recall 
that f3 > 0). 

If Mi >r (s) has simple zeroes (so s 7^ Si) then, similarly, the flux out of 
the left/right sides is given by: 
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a e b 

Figure 16: The Lyapunov dimension variation. 

a) Samples of e values along the bifurcation curve e = sj r (uj, 7, /3, 6; £i r ) with 
fj = 0, i = 0.05 and 4. =0. 

b) Zoom in on e values in the interval [0.2,0.25]. 



= r i M l , r (s)ds= (28) 



so 



1 — ( — j + — arcsm( 



— (~iFi t k 1 — I I -f — arcsmi 



In region II, A° ut = 0, hence obviously |A°"*| > |A°"*| in this region, and it 
can be shown that |A°"*| > \A° ut \ in region III as well (i.e. for 6 < ^- L ). 

One might expect that the ratio between the fluxes to the right and left 
regions determines the ratio between the sizes of the basins of attraction. 
However, in regions II and III near the boarder line between the regions, this 
picture may change dramatically; There are cases for which all the initial 
conditions which are numerically integrated are attracted to the right side. 

Numerical calculations of the basins of attraction of the left /right attrac- 
tors suggest a more detailed description: 

1. No intersections on both sides - region I: 

To order e, by equation (27) the left and right influxes equal, 
but high order terms alter these results. Indeed, numerically it is found 
that the right basin is larger than the left one. Moreover, it's area seems 
to grow monotonically as 6 decreases to it's threshold value, 37-F//4. 
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Figure 17: The relative size of the right and left basins of attractions. 

— : the right basin size, : the left basin size, *: sampled values of e. 

The secondary homoclinic bifurcation values e(£ cc = 1) are denoted by 'x', 
and the secondary homoclinic bifurcation values e(£ c d = 1) are denoted by 



V. 



2. Intersections on the left side or on both sides - regions II and III: 

As 6 is further decreased, the area of the right basin continues to grow. 
However, the growth in the basin's area seems to be discontinuous. 
This phenomenon is associated with the "boundary metamorphosis" 
[15] of the subharmonics. For some parameter values in region III a 
two sided periodic orbit or SA may exist, hence in these cases, some 
(or most) of the initial conditions are not attracted neither to a left nor 
to a right attractor. 

There seems to be a correlation between the discontinuities and the occur- 
rence of homoclinic tangencies. See for example figure 17. In this figure the 
percentage of the sampled phase space area which is swept to the left/right 
sides, up to an approximated error of ±0.03, is presented as a function of e 
(all the other parameters are fixed), were the values of e which correspond 
to secondary homoclinic bifurcation curves are specified. The results of this 
figure are somewhat puzzling; One would expect that left/right attractors 
are composed, for £ > 0, from the attracting resonances. Hence, when the 
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unstable manifold intersects through the resonance region, it depletes the 
resonance, thus decreasing the relative area of the basin of attraction of the 
corresponding side. However, in figure 17, along with results that confirm 
this scenario (see for example the jump near e = 0.06), we observe quite 
the contrary results (see for example the huge jump near e = 0.5). Possibly 
other, undetected bifurcation is responsible to these results. 

In [21], numerical results regarding the relation between homoclinic and 
other bifurcation curves and the basins of attraction of systems with a cubic 
like potential well are presented. There, it has been suggested that the bi- 
furcation curve, corresponding to what we call here a secondary homoclinic 
bifurcation curve with a structural index i = 0, is of great significance, since 
closely beneath it they numerically observed a chaotic escape (i.e. a destruc- 
tion of the basin of attraction of the SA). Their chaotic escape corresponds, 
in a case of a closed 11 system, to a decrease in the basin of attraction of one 
side and an increase in the basin of attraction of the other side. Thus, the 
current results are in agreement with the results obtained in [21] for open 
systems. Possibly, the critical curves (near which SA appear or loose stabil- 
ity) that they have observed numerically, correspond to the curves discussed 
above: the homoclinic bifurcation curve above which the lobe El\ intersects 
the lobe Dl at at-least five points, and the curve at which the unstable 
manifold intersects the resonance region. 

4 Discussion and conclusions 

The qualitative differences between the flows under symmetric and asym- 
metric forcing loom when primary homoclinic intersections/tangencies occur 
only on one side of the saddle fixed point (region II + its neighborhood). 
This region may be of significant size even for very small asymmetry values 
( /? <C 1) if the forcing frequency is appropriately set. It is of negligible size 
in the adiabatic limit, hence, to the best of our knowledge, was not observed 
in previous works which have considered asymmetric potentials with adia- 
batic forcing. In this fat region II, the system may posses one one-sided-SA 
(strange attractors), two one-sided-SA, or one two-sided SA. In the former 
case the attractor may be situated on either side of the fixed point; however, 

11 A system is called closed if some forward iteration of the Poincare map of a segment 
of the unstable manifold which has left the left/right region returns to it. See [20] for more 
precise definition. 
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the nature of the basin of attraction of the left and right SA seems to be 
different. 

We find that the relative size of the basin of attraction to the left /right 
attractors is usually not sensitive to its strangeness (i.e. the size of the basin 
does not change significantly when a SA is destroyed/created). In general, 
the basin of attraction of the right attractor is always a bit larger than that 
of the left attractor, where in most cases the difference between the fluxes 
to the right and left regions determines the ratio between the sizes of the 
basins of attraction. This occurs in a continuous and natural way in region I, 
However, near the border between regions II and III, where one sided SA may 
appear (on either side of it), this picture may change dramatically; There are 
cases for which all the initial conditions which are numerically integrated are 
attracted to the right side. 

The robust, observable (hence physically significant) SA appear near pri- 
mary homoclinic tangencies and near secondary homoclinic tangencies with 
small structural index. A key perturbational tool for finding the latter is 
the SMF (secondary Melnikov function [33]). Generally, it is found that the 
SMF supplies excellent analytical prediction to the occurrence of secondary 
homoclinic tangencies even for relatively large values of e. However, it fails 
near the £ = homoclinic bifurcations, exactly in the region were robust SA 
exist. Thus, only lower bounds and approximate curves for the regions were 
SA are observed are found. We suspect that both phenomena (the failure of 
the SMF and the appearance of SA) are associated with the involvement of 
the stable and unstable manifolds in a I : I resonance. Thus we derived a 
simple lower bound for e above which the manifolds enter the I : I resonance 
was derived (eq. (25)). The study of the relation between the resonance, the 
manifolds and the strange attractor, and the construction of a more accurate 
approximation to the homoclinic bifurcation curves near a I : I resonance 
are left for future work. 

The structure of the SA vary with the parameters; as the ratio ^ increases 
the values of the positive Lyapunov exponent slightly decreases and the Lya- 
punov (fractal) dimension decreases significantly. Surprisingly, we find that 
fixing this ratio and varying the other parameters in one of the "windows" for 
which SA exist, the structure, Lyapunov exponent and Lyapunov dimension 
of the SA hardly change. Such a variation does change, in particular, the 
dissipation ( contracting) rate per Poincare map. 
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A Finding secondary homoclinic bifurcation 
curves. 

In this appendix some technical aspects regarding the method of solution of 
( 17) - (21) are described. 

From (17) - (21) we construct the equations: 



m: 



-M c (t j) + Jcd T, 



(29) 



{ P -\t\Jt ) - t ) 

M c {t ) 
P-\2(t\ cd (t ) - t )) 

M c {t ) 



c = d 



(30) 



c ^ d 



where c, d £ {/, r}; i = 1, 2; j cd £ N, and from equation (7) for the Melnikov 
function of the AFDO, one gets: 



m: 



M^' 2 (x) 



Notice that Mj 1,1 (x) are undefined for 



1 . x U 

— arcsml — — — — H — — — ) , 

TV 1 X 4:6 

arcsml — — — — H — — — 

to to K lF d (u>,P) 3 7 F d {ioJ)> 



+ 
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> 1 , hence there 



are some values of to in [0, T) for which e l cd {to) from (30) are undefined. From 

(20) and (23) we get that for to £ [0,T): j cd = £ cd + s(t ), where s(t ) = 
for t £ [0,^-), and s(t ) = 1, for to £ [|-,T). The above equations ((29) 
and (30)) to be dependent on the perturbation parameters: f3 (asymmetry), 
6 (dissipation) 7 (the amplitude of the forcing) and u (the frequency of the 
for forcing). In addition, for \H\ = \eM c {to)\ <C 1 we get (see equations (6), 

(21) , (22) and (29)): 



dhf(t ,e) 
dt 



M' c (to) + M' d (M^(-M c (to)))(l 



M' c (t ) 

M c {t ) 



(31) 
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c, d £ {/, r}, z = 1,2. 
Hence, for the AFDO, equation (29) becomes: 

*L *o = - 27r + arcsin — sin w f + — — — - ] 32 

2tt 

+ (icd + s(t )) — , 

2 1 F c (lo } /3) 8<5 

*icd(*o) = -[vr-arcsm -— — sin cut + — — -~)\ 33 

l< d (LO,8) Sjl< d (LLi,8) 

+ (icd + S(t )) — , 

f 0, t o G[0,-) 

S(* ) = < vr^TT 5 C,G?£{/,r}, 

1, t G [-,— ) 

and equation (30) becomes: 

f -16exp(* -<U*o)) 

4d(*o) ~ 



7 F c (cu) sin(^o) - f ' 
16exp(t -t\ cd {t )) 



7 F c (u;)sin(u;*o)-f ' C ^ d ' 
i = 1,2. (34) 

These approximations are valid for sufficiently small H's of the period func- 
tion (see equation (6)), namely for: 

4<5 

4 c (t )[ 7 FcHsin(^o)-y] ^ 0- (35) 

(36) 

4(5 

^(M^cHsin^-y] ^ + 



(37) 
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And, for such sufficiently small values of H } equation (31) is: 

— — sin (ut ) + cos (ut\ cd (t [tan (ut ) 

to sin(cuto) 



(38) 



sin(cuto) 



3 7 F c (coJ) 



where F C} are as in (9), (10) of section 2.1, for c, d £ 0,^}, and £ l cd (t ) is 
calculated here with the use of the approximated value of the period function 
of AFDO, P(H), from (6). 

Remark 1 Since the exact inverse function of the period function of the 
AFDO, P~ 1 (x), cannot be found analytically, (t 0} e^(^o))j * = 1 , 2 are found 
by solving the equations: 

P(eM c (t )) = r^to) (39) 
dP(eM c (t )) drU{to) 



dt dt 



(40) 



wh 



ere: 



r cd {t ) - | 2{t[Jto) _ to); Mc[to) > ^ ^ , 141 j 

c,de{l,r}, i G {1,2}, t o e[0,T), T = — , 

u 

by a Newton method, combined with a linear prolongation scheme. As the 
initial guesses for the Newton method we use the approximated values for 
t and e l cd {t l ). See [20] for details on how this approximated values may be 
obtained. 

When approximation to leading order in H for P(H) (as in equation (6)) 
is used to calculate £ l cd (t Q ), one gets that to leading order in f3, |e^(^o)l = 
l e cc(to)l for 6 = 0. Actually, more accurate approximations (as using higher 
terms and Newton method) show that they are different. 
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